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^13 Summary 

^4 We consider variable selection in high-dimensional linear models where the number of covari- 

^ 15 ates greatly exceeds the sample size. We introduce the new concept of partial faithfulness and use 

I— ,16 it to infer associations between the covariates and the response. Under partial faithfulness, we 

17 develop a simplified version of the PC algorithm (Spirtes et al., 2000), the PC-simple algorithm, 

^ which is computationally feasible even with thousands of covariates and provides consistent vari- 
able selection under conditions on the random design matrix that are of a different nature than 

^ 20 coherence conditions for penalty-based approaches like the Lasso. Simulations and application 

t^Jl to real data show that our method is competitive compared to penalty -based approaches. We 

22 provide an efficient implementation of the algorithm in the R-package pcalg. 

m23 

^ 24 Some key words: Directed acyclic graph; Elastic net; Graphical modeling; Lasso, Regression 

cn27 

^^28 1- Introduction 

O 29 Variable selection in high-dimensional models has recently attracted a lot of attention. A par- 

^ 30 ticular stream of research has focused on penalty-based estimators whose computation is feasible 

. . 31 and provably correct (Meinshausen & Biihlmann, 2006; Zou, 2006; Zhao & Yu, 2006; Candes 

.^32 & Tao, 2007; van de Geer, 2008; Zhang & Huang, 2008; Wainwright, 2009; Meinshausen & 

^ 33 Yu, 2009; Huang et al., 2008; Bickel et al., 2009; Wasserman & Roeder, 2009; Candes & Plan, 

^ $4 2009). Another important approach for estimation in high-dimensional settings, including vari- 

- - 35 able selection, has been developed within the Bayesian paradigm, see for example George & Mc- 

36 Culloch (1993, 1997); Brown et al. (1999, 2002); Nott & Kohn (2005); Park & Casella (2008). 

37 These methods rely on Markov chain Monte Carlo techniques which are typically very expensive 

38 for truly high-dimensional problems. 

39 In this paper, we propose a method for variable selection in linear models which is funda- 

40 mentally different from penalty-based schemes. Reasons to look at such an approach include: 

41 (i) From a practical perspective, it is valuable to have a very different method in the tool-kit 

42 for high-dimensional data analysis, raising the confidence for relevance of variables if they are 

43 selected by more than a single method, (ii) From a methodological and theoretical perspective, 

44 we introduce the new framework of partial faithfulness. Partial faithfulness is related to, and typ- 

45 ically weaker than, the concept of linear faithfulness used in graphical models, hence the name 

46 partial faithfulness. We prove that partial faithfulness arises naturally in the context of linear 

47 models if we make a simple assumption on the structure of the regression coefficients to exclude 

48 adversarial cases, see assumption (C2) and Theorem 1 . 
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49 Partial faithfulness can be exploited to construct an efficient hierarchical testing algorithm, 

50 called the PC-simple algorithm, which is a simplification of the PC algorithm (Spirtes et al., 

51 2000) for estimating directed acychc graphs. The letters PC stand for the first names of Peter 

52 Spirtes and Clarke Glymour, the inventors of the PC algorithm. We prove consistency of the 

53 PC-simple algorithm for variable selection in high-dimensional partially faithful Unear models 

54 under assumptions on the design matrix that are very different from coherence assumptions for 

55 penalty-based methods. The PC-simple algorithm can also be viewed as a generalization of cor- 

56 relation screening or sure independence screening (Fan & Lv, 2008). Thus, as a special case, we 

57 obtain consistency for correlation screening under different assumptions and reasoning than Fan 

58 & Lv (2008). We illustrate the PC-simple algorithm, using our implementation in the R-package 

59 pcalg, on high-dimensional simulated examples and a real data set on riboflavin production by 

60 the bacterium Bacillus subtiUs. 
61 

62 

63 2. Model and notation 

64 Let X = (X(^), . . . , G be a vector of covariates with E{X) = fix and cov{X) = 

65 T,x- Let e G M with E{e) = and var(e) = cr^ > 0, such that e is uncorrected with 

66 X^^^ , . . . ,X^\ Let y G M be defined by the following random design linear model: 

68 
69 
70 



y = 5 + ^/3^.xW+e, (1) 

3=1 



71 for some parameters 5 G R and /? = (/?i, . . . , /3p)^ G W. We assume implicitly that E{Y'^) < 

72 oo and E{{X(^^f} < oo for j = 1, . . . ,p. 

73 We consider models in which some, or most, of the Pjs are equal to zero. Our goal is to identify 

74 tiie active set 

76 A = {j = l,...,p; P,^0} 

77 based on a sample of independent observations {Xi,Yi), . . . , (X„, ¥„) which are distributed as 

78 (X, y). We denote the effective dimension of the model, that is, the number of nonzero PjS, by 

79 peff =1^1. We define the following additional conditions: 
80 

81 (CI) Tix is strictly positive definite. 
82 

83 (C2) The regression coefficients satisfy 
84 
85 

86 where /(•) denotes a density on a subset of MP^^ of an absolutely continuous distribution 

87 with respect to Lebesgue measure. 
88 

§9 Assumption (CI) is a condition on the random design matrix. It is needed for iden- 

9Q tifiability of the regression parameters from the joint distribution of {X,Y), since (3 = 

91 I]^^cov(y, . . . , cov(y, . Assumption (C2) says that the non-zero regression co- 

92 efficients are realizations from an absolutely continuous distribution with respect to Lebesgue 

93 measure. Once the 13 js are realized, we fix them such that they can be considered as deterministic 

94 in the linear model (1). This framework is loosely related to a Bayesian formulation treating the 

95 PjS as independent and identically distributed random variables from a prior distribution which 

96 is a mixture of a point mass at zero for PjS with j ^ A and a density with respect to Lebesgue 



{Pj; jeA}^f{b)db, 
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97 measure for PjS with j G A. Assumption (C2) is mild in the following sense: the zero coeffi- 

98 cients can arise in an arbitrary way and only the non-zero coefficients are restricted to exclude 

99 adversarial cases. Interestingly, Candes & Plan (2009) also make an assumption on the regres- 

100 sion coefficients using the concept of random sampUng in their generic S-sparse model, but other 

101 than that, there are no immediate deeper connections between their setting and ours. Theorem 1 

102 shows that assumptions (CI) and (C2) imply partial faithfulness, and partial faithfulness is used 

103 to obtain the main results in Theorems 3, 4 and 5. Assumption (C2), however, is not a necessary 

104 condition for these results. 

105 We use the following notation. For a set 5 C {1, . . . \S\ denotes its cardinality, 5*^ is 

106 its complement in {1, . . . ,p}, and X'-^^ = {X^^); j G S}. Moreover, Z^^) | W) and 

107 parcov(Z(^), Z^^) | W) denote the population partial correlation and the population partial co- 

108 variance between two variables Z^^^ and Z^^^ given a collection of variables W. 
109 

110 

111 3. Linear FAITHFULNESS AND PARTIAL FAITHFULNESS 

112 3-1. Partial faithfulness 

113 We now define partial faithfulness, the concept that will allow us to identify the active set A 

114 using a simplified version of the PC algorithm. 
115 

llf^ Definition 1. Let X €zMP be a random vector, and let Y be a random variable, 

llj The distribution of {X,Y) is said to be partially faithful if the following holds for every 

118 J G {1, . . . ,p}: ifp{Y, XO) I = Ofor some S C {j}^ then p(Y, X^^) \ Xdj}"")) = 0. 

120 
121 

122 p{Y, X(^) I X^'^)) = for some <S C {j}^ implies pj = 0. (2) 
123 

124 Theorem 1 . Consider the linear model (1) satisfying assumptions (CI) and ( C2). Then par- 

125 tial faithfulness holds almost surely with respect to the distribution generating the non-zero re- 

126 gression coefficients. 

127 A proof is given in the Appendix. This is in the same spirit as a result by Spirtes et al. (2000, 

128 Th. 3.2) for graphical models, saying that non-faithful distributions for directed acyclic graphs 

129 have Lebesgue measure zero, but we are considering here the typically weaker notion of partial 

130 faithfulness. A direct consequence of partial faithfulness is as follows: 

131 

132 Corollary 1. Consider the linear model (1) satisfying the partial faithfulness condition. 

133 Then the following holds for every j E {1, . . . ,p}: 

J^^ p{y, X^^^ I X^^^) + Ofor all S C {j}^ if and only if/3j + 0. 

136 A simple proof is given in the Appendix. Corollary 1 shows that, under partial faithfulness, 

137 variables in the active set A have a strong interpretation in the sense that all corresponding 

138 partial correlations are different from zero when conditioning on any subset <S C {j}*^. 
139 

140 3-2. Relationship between linear faithfulness and partial faithfulness 

141 To clarify the meaning of partial faithfulness, this section discusses the relationship between 

142 partial faithfulness and the concept of linear faithfulness used in graphical models. This is the 

143 only section that uses concepts from graphical modeling, and it is not required to understand the 

144 remainder of the paper. 



For the linear model (1) with assumption (CI), /3j = Oif andonly if | x({^>^)) = 0. 

Hence, such a model satisfies the partial faithfulness assumption if for every j E {1, . . . ,p}: 
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145 We first recall the definition of linear faithfulness. The distribution of a collection of random 

146 variables can be depicted by a directed acyclic graph G in which each vertex 

147 represents a variable, and the directed edges between the vertices encode conditional depen- 

148 dence relationships. The distribution of{Z^^\... , ) is said to be linearly faithful to G if the 

149 following holds for alH 7^ j G {1, . . . , and 5 C {1, . . . , q} \ {i, j}: 

j^^ Z^^ and Z^^ are d-separated by Z^^^ in G if and only if p{Z^^ , Z^^ \ Z^^^) = 0, 

152 see, e.g., Spirtes et al. (2000, page 47). In other words, linear faithfulness to G means that all and 

153 only all zero partial correlations among the variables can be read off from G using d-separation, 

154 a graphical separation criterion explained in detail in Spirtes et al. (2000). 

155 Partial faithfulness is related to a weaker version of linear faithfulness. We say that the distri- 

156 bution of {X, Y), where X G is a random vector and y G M is a random variable, is linearly 

157 r-faithful to G if the following holds for all j G {1, . . . , p} and 5 C {jjC; 
158 

and Y are d-separated by X^^^iaG if and only if p{X^^\Y \ X^^^) = Q. (3) 

160 Thus, linear F-faithfulness to G means that all and only all zero partial correlations between Y 

161 and the X^^\ can be read off from G using d-separation, but it does not require that all and only 

162 all zero partial correlations among the X'^^h can be read off using d-separation. 

163 We now consider the relationship between linear faithfulness, linear F-faithfulness and par- 

164 tial faithfulness. Linear faithfulness and linear y-faithfulness are graphical concepts, in that they 

165 link a distribution to a directed acyclic graph, while partial faithfulness is not a graphical con- 

166 cept. From the definition of linear faithfulness and linear y-faithfulness, it is clear that linear 

167 faithfulness impUes Unear F-faithfulness. The following theorem relates Unear F-faithfulness to 

168 partial faithfulness: 
169 

Theorem 2. Assume that the distribution of {X,Y) is linearly Y-faithful to a directed 
acyclic graph in which Y is childless. Then partial faithfulness holds. 

172 A proof is given in the Appendix. A distribution is typically linearly F-faithful to several 

173 directed acyclic graphs. Theorem 2 applies if Y is childless in at least one of these graphs. 
174 

We illustrate Theorem 2 by three examples. Example 1 shows a distribution where partial 

^rj^ faithfulness does not hold. In this case. Theorem 2 does not apply, because the distribution of 
{X, Y) is not linearly y-faithful to any directed acyclic graph in which Y is childless. Examples 
2 and 3 show distributions where partial faithfulness does hold. In Example 2, the distribution of 

^rjg {X, Y) is Hnearly y-faithful to a directed acyclic graph in which Y is childless, and hence partial 
faithfulness follows from Theorem 2. In Example 3, the distribution of {X, Y) is not linearly Y- 

J faithful to any directed acyclic graph in which Y is childless, showing that this is not a necessary 

^ condition for partial faithfulness. 

183 Example 1. Consider the following Gaussian linear model: 

j^^ = ei, = + £2, y = XW - + e, (4) 

186 where ei, 82 and e are independent standard Normal random variables. This model can be rep- 

187 resented by the linear model (1) with /?i = 1 and /?2 = — 1. Furthermore, the distribution of 

188 {X, Y) = (X(^) , X(^) , y) factorizes according to the graph in Figure 1(a). 

189 The distribution of {X, Y) is not partially faithful, since p{Y, | 0) = but p{Y, X^^^ \ 

190 / 0. Theorem 2 does not apply, because the disttibution of {X,Y) is not hnearly Y- 

191 faithful to any directed acyclic graph in which Y is childless. For instance, the distribution of 

192 (X, y) is not hnearly y-faithful to the graph in Figure 1(a), since /9(X(^) , y | 0) = but X(^) 
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193 and Y are not d-separated by the empty set. The zero correlation between X^^^ and Y occurs 

194 because X^^^ = ei drops out of the equation for Y due to a parameter cancellation that is simi- 

195 lar to equation (Al) in the proof of Theorem 1:1" = X^^) - X^^) + ^ = _ (e^ + ^2) + e = 

196 + £. The distribution of {X, Y) is linearly faithful, and hence also linearly Y-faithful, to the 

197 graph X^^^ X(^) <— Y, but this graph is not allowed in Theorem 2 because Y has a child. 

198 Such failure of partial faithfulness can also be caused by hidden variables. To see this, consider 

199 the following Gaussian hnear model: 
200 

201 = £1, = £3, = + +e2,Y = + e, 
202 

203 where ei,e2,£3 and e are independent standard Normal random variables. The distribution of 

204 , , , Y) factorizes according to the DAG ^ ^ ^ Y, and is lin- 

205 early faithful to this DAG. Hence, the distribution of (X^^) , X^^) , X^^) , Y) is partially faithful by 

206 Theorem 2. If, however, variable X^^^ is hidden, so that we only observe y), then 

207 the distribution of (X^^^ , X(^) , Y) has exactly the same conditional independence relationships 

208 as the distribution arising from (4). Hence, the distribution of {X^^\X^'^\Y) is not partially 

209 faithful. 
210 

21 1 Example 2. Consider the following Gaussian hnear model: 

212 

213 = £1, = xw + £2, = x(^) + £3, ^('^ = x^^^ - + £4, y = + £, 

214 

215 where £1, . . . , £4 and £ are independent standard Normal random variables. This model can be 

216 represented by the linear model (1) with /3i = /33 = /34 = and /32 = 1. Furthermore, the distri- 

217 bution of {X, Y) = {X^^\ . . . , X^'^\ Y) factorizes according to the graph in Figure 1(b). 

218 The distribution of {X, Y) is partially faithful, since p{Y, X^) | x({-'>'')) ^ only for j = 2, 

219 and p(y, X^^) | X^"^)) 7^ for any S C {1,3,4}. In this example, partial faithfulness follows 

220 from Theorem 2, since the distribution of {X, Y) is linearly Y-faithful to the graph in Figure 1(b) 

221 and Y is childless in this graph. The distribution of {X, Y) is not linearly faithful to the graph 

222 in Figure 1(b), since cor{X^^\ X^^^) = but X^^^ and X^^^ are not d-separated by the empty 

223 set. Moreover, there does not exist any other directed acyclic graph to which the distribution 

224 of {X, Y) is linearly faithful. Hence, this example also illustrates that linear Y-faithfulness is 

225 strictly weaker than linear faithfulness. 
226 

227 Example 3. Consider the following Gaussian hnear model: 
228 

229 X(i) = £1, = + £2, = + £3, y = - + £, 
230 

231 where £1, £2, £3 and e are independent standard Normal random variables. This model can be 

232 represented by the linear model (1) with /3i = 0, /32 = 1 and = —1- Furthermore, the distri- 

233 bution of {X, Y) = {X^^'^ , X^'^^ , X^^^ , Y) factorizes according to the graph in Figure 1(c). 

234 The distribution of {X, Y) is partially faithful, since p{Y, X^^) | X^^^^''^) / for j € {2, 3}, 

235 p(Y, I X(^^) ^ for any S C {1, 3}, and p{Y, X(3) | x('5)) / for any 5 C {1, 2}. How- 

236 ever, in this case partial faithfulness does not follow from Theorem 2, since the distribution of 

237 (X, Y) is not linearly F-faithful to the graph in Figure 1(c), since cor (X(i),F) = but X(y 

238 and Y are not d-separated by the empty set. Moreover, there does not exist any other directed 

239 acychc graph to which the distribution of (X, Y) is linearly y-faithful. 
240 
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(a) Example 1 
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(b) Example 2 



X{3) ^ Y 

(c) Example 3 



Fig. 1. Graphical representation of the models used in Ex- 
amples 1-3. 

4. The PC-simple algorithm 
41. Population version of the PC-simple algorithm 
We now explore how partial faithfulness can be used for variable selection. In order to show 
the key ideas of the algorithm, we first assume that the population partial correlations are known. 
In Section 4-2 we consider the more realistic situation where partial correlations are estimated 
from data. 

First, using 5 = in expression (2) yields that Pj = if cor{Y, X^j^) = for some j G 
{1, . . . ,p}. This shows that the active set A cannot contain any j for which cor(y,X(j)) = 0. 
Hence, we can screen all marginal correlations between pairs (Y, X^^^), j = 1,. . . ,p, and build 
a first set of candidate active variables 



A^^] = {j = l,...,p; cor(y,x(^))7^0}. 



(5) 



We call this the step^ active set or the correlation screening active set, and we know by partial 
faithfulness that 



(6) 



Such correlation screening may reduce the dimensionality of the problem by a substantial 
amount, and due to (6), we could use other variable selection methods on the reduced set of 
variables A^^\ 

Furthermore, for each j G A^^^ expression (2) yields that 

p{Y, X(^) I X^''^) = for some k G A^^^ \ {j} impUes /3j = 0. (7) 

That is, for checking whether the jth covariate remains in the model, we can additionally screen 
all partial correlations of order one. We only consider partial correlations given variables in the 
stepx active set A^^^ This is similar to what is done in the PC algorithm, and yields an important 
computational reduction while still allowing us to eventually identify the true active set A. Thus, 
screening partial correlations of order one using (7) leads to a smaller active set 

^[2] = {j e ^[1]; p(y,x(^) I for all k e \ {j}} c ^W. 

This new step2 active set A^'^^ further reduces the dimensionality of the candidate active set, 
and because of (7) we still have that A^^^ ^ A. We can continue screening higher-order partial 
correlations, resulting in a nested sequence of step^ active sets 

^[1] 2 ^[2] D...D D---DA. (8) 

A stcp„ active set A^"^^ could be used as dimension reduction together with any favored vari- 
able selection method in the reduced linear model with covariates corresponding to indices in 
. Alternatively, we can continue the algorithm until the candidate active set does not change 
anymore. This leads to the PC-simple algorithm, shown in pseudocode in Algorithm 1 . 
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289 Algorithm 1 The population version of the PC-simple algorithm. 

•^^^ 1: Set m = 1. Do correlation screening, and build the stepi active set 

291 ^[1] = {j = 1, . . . ,p; cor(y, XO)) 7^ 0} as in (5). 

^92 2: repeat 

3. m = m + 1. Construct the step™ active set: 

294 

295 ^[""1 = {j € ; piY, X^^^ \ X^'^)) / 

296 for all <S C ^1™"!] \ {j} with \S\ = m - 1}. 

297 ~ \\JS \ \ / 

298 4: until \A^'^'^\ < m. 

299 



300 
301 



319 
320 
321 



325 



The value m that is reached in Algorithm 1 is called mjcach: 



302 rwreach = min{m; l^l^lj < m}. (9) 

303 

The following theorem shows correctness of the population version of the PC-simple algorithm. 

305 Theorem 3. For the linear model (1) satisfying (CI) and partial faithfulness, the population 

306 version of the PC-simple algorithm identifies the true underlying active set, i.e., ["breach] = _4 = 

307 {j = l,...,p; Pjj^O}. 

A proof is given in the Appendix. Theorem 3 shows that partial faithfulness, which is often 
weaker than linear faithfulness, is sufficient to guarantee correctness of the population PC-simple 
algorithm. The PC-simple algorithm is similar to the PC algorithm (Spirtes et al., 2000, Section 
5-4-2), but there are two important differences. First, the PC algorithm considers all ordered 
pairs of variables in (X^^^ , . . . , X^^ , Y), while we only consider ordered pairs (y,X(^)), j e 
{1, . . . ,p}, since we are only interested in associations between Y and X^^\ Second, the PC 
algorithm considers conditioning sets in the neighborhoods of both Y and X^^\ while we only 
consider conditioning sets in the neighborhood of Y. 



308 
309 
310 
311 
312 
313 
314 
315 
316 

317 4-2. Sample version of the PC-simple algorithm 

° For finite samples, the partial correlations must be estimated. We use the following shorthand 

notation: 

PiY, j\S)= PiY, X(^-) I ) , PiY, j\S)= PiY, X(-?) I ) , 

322 PiiJ I 3) = PiX^'Kx^^) I X(^)), pii,j I S) = | X^^)), 

323 where the hat-versions denote sample partial correlations. These can be calculated recursively, 

324 since for any G (S we have 

piY,j I S \ {k}) -piY,k\S\ {k})pij, k\S\ {k}) 



fll ^^^'^ ' [{1 -p{Y.k\S\ {fc})2}{l - pij, k\S\ {fe})2}] V2 • 

328 In order to test whether a partial correlation is zero, we apply Fisher's .Z-transform 

330 Z(y,;|S) = ilog{l±4Mi|l|. (10) 

331 2 yi- piY,j\S)] 

332 Classical decision theory in the Gaussian case yields the following rule. Reject the nuU- 

333 hypothesis H^iY,] \ S) : piY,j | 5) = against the two-sided alternative HAiY,j \ S) : 

334 piY,j I 5) 7^ if (ra - |5| - 3y/'^\ZiY,j \ S)\ > $"^(1 - a/2), where a is the significance 

335 level and $(•) is the standard Normal cumulative distribution function. Even in the absence of 

336 Gaussianity, the rule above is a reasonable thresholding operation. 
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337 The sample version of the PC-simple algorithm is obtained by replacing the statements about 

338 p{Y, XO) I ) ^ in Algorithm 1 by 



352 
353 
354 



368 
369 



_ \s\ - 3)V2|z(y, j I 5)1 > - a/2). 



339 
340 

^■^1 The resulting estimated set of variables is denoted by A{a) = A"^''^'^''^{a), where mreach is the 

^42 estimated version of the quantity in (9). The only tuning parameter a of the PC-simple algorithm 

is the significance level for testing the partial correlations. 

The PC-simple algorithm is very different from a greedy scheme, since it screens many corre- 
lations or partial correlations at once and may delete many variables at once. Furthermore, it is a 
more sophisticated pursuit of variable screening than the marginal correlation approach in Fan & 
Lv (2008) or the low-order partial correlation method in Wille & Buhlmann (2006). Castelo & 
Roverato (2006) extended the latter and considered a limited-order partial correlation approach. 
However, their method does not exploit the clever trick of the PC-simple algorithm that it is suf- 
ficient to consider only conditioning sets S which have survived in the previous step^_i active 
set ^["^"^1. Therefore, the algorithm of Castelo & Roverato (2006) is often infeasible and must 
be approximated by a Monte Carlo approach. 

Since the PC-simple algorithm is a simplified version of the PC algorithm, its computational 
complexity is bounded above by that of the PC algorithm. This is difficult to evaluate exactly, but 
a crude bound is 0{p^'^^), see Kalisch & Buhlmann (2007, formula (4)). Section 6 shows that 
we can easily use the PC-simple algorithm in problems with thousands of covariates. 

358 

359 5. Asymptotic results in high dimensions 

5 1. Consistency of the PC-simple algorithm 
We now show that the PC-simple algorithm is consistent for variable selection, even if p is 
much larger than n. We consider the linear model in (1). To capture high-dimensional behavior, 
364 dimension grow as a function of sample size and thus, p = Pn and also the distribution 

of {X, y), the regression coefficients (3j = Pj^n^ and the active set A = An with peff = pefF^ = 
366 l^^l change with n. Our assumptions are as follows: 

(Dl) The distribution P„ of (X, Y) is multivariate Normal and satisfies (CI) and the partial faith- 
fulness condition for all n. 
(D2) The dimension p„ = 0(n") for some < a < oo. 
(D3) The cardinality of the active set pefF„ = \ An\ = |{i = Ij • • • ,Pn; Pi,n 7^ 0}| satisfies: 
3'71 pcff„ = 0{n^-^) for some < 5 < 1. 

372 (D4) The partial correlations p„(y, j\S)= p{Y, | X('5)) satisfy: 

374 inf ||p„(y, j \S)\-j = l,... ,pn, S C {jf, \S\ < pefF„ with pn{Y,j \S)^o}>Cn, 

375 

376 where = 0{n'^) for some < d < h/2, and h is as in (D3). 

377 (D5) The partial correlations pniY,j \ S) and pn(i,3 \ S) = | X^-^)) satisfy: 
378 
379 
380 
381 
382 

383 Assumption (Dl) is made to simplify asymptotic calculations, and it is not needed in the popu- 

384 lation case. Unfortunately, it is virtually impossible to check assumptions (D1)-(D5) in practice. 



(i) sup |p„(y,j |5)| <M< 1, 

n,j,SC{j}C,\S]<peS„ 

(m) sup \pn{i,j \S)\<M <1. 

n,i^j,SC{i,j}C,\S\<pef£^ 
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385 with the exception of (D2). However, this is common to assumptions for high-dimensional vari- 

386 able selection, such as the neighborhood stability condition (Meinshausen & Buhlmann, 2006), 

387 the irrepresentable condition (Zhao & Yu, 2006), or the restrictive eigenvalue assumption (Bickel 

388 et al., 2009). A more detailed discussion of assumptions (D1)-(D5) is given in Section 5-2. 

389 Letting An{c() denote the estimated set of variables from the PC-simple algorithm in Section 

390 4-2 with significance level a, we obtain the following consistency result: 
391 

Theorem 4. Consider the linear model (1) and assume (D1)-(D5). Then there exists a se- 
quence — >^ (n ^ oo) and a constant C > such that the PC-simple algorithm satisfies 

394 pr{An{an) = An} = l- Oie^-Cn^-^'^)} ^ I {n ^ oo), 
395 

where a is as in (D4). 

397 A proof is given in the Appendix. The value a„, although being the significance level of a single 

398 test, is a tuning parameter which allows to control type I and 11 errors over the many tests which 

399 are pursued in the PC-simple algorithm. A possible choice yielding consistency is = 2{1 — 

400 $(n^/^c„/2)}. This choice depends on the unknown lower bound of the partial correlations in 

401 (D4). 
402 

403 5-2. Discussion of the conditions of Theorem 4 

404 There is much recent work on high-dimensional and computationally tractable variable se- 

405 lection, most of it considering versions of the Lasso (Tibshirani, 1996) or the Dantzig selector 

406 (Candes & Tao, 2007). Neither of these methods exploit partial faithfulness. Hence, it is interest- 

407 ing to discuss our conditions with a view towards these other established results. 

408 For the Lasso, Meinshausen & Buhlmann (2006) proved that a so-called neighborhood sta- 

409 bility condition is sufficient and almost necessary for consistent variable selection, where the 

410 word almost refers to the fact that a strict inequality with the relation < appears in the suffi- 

41 1 cient condition whereas for necessity, there is a < relation. Zou (2006) and Zhao & Yu (2006) 

412 gave a different, but equivalent condition. In the latter work, it is called the irrepresentable con- 

413 dition. The adaptive Lasso (Zou, 2006) or other two-stage Lasso and thresholding procedures 

414 (Meinshausen & Yu, 2009) yield consistent variable selection under weaker conditions than the 

415 neighborhood stability or irrepresentable condition, see also Example 4 below. Such two-stage 

416 procedures rely on bounds for ||/3 — P\\q {q = 1, 2) whose convergence rate to zero is guaranteed 

417 under possibly weaker restricted eigenvalue assumptions on the design (Bickel et al., 2009) than 

418 what is required by the irrepresentable or neighborhood stability condition. All these different 

419 assumptions are not directly comparable with our conditions (D1)-(D5). 

420 Assumption (D2) allows for an arbitrary polynomial growth of dimension as a function of 

421 sample size, while (D3) is a sparseness assumption in terms of the number of effective variables. 

422 Both (D2) and (D3) are fairly standard assumptions in high-dimensional asymptotics. More crit- 

423 ical are the partial faithfulness requirement in (Dl), and the conditions on the partial correlations 

424 in (D4) and (D5). 

425 We interpret these conditions with respect to the design X and the conditional distribution of 

426 Y given X. Regarding the random design, we assume (CI) and (D5,(ii)). Requiring (CI) is rather 

427 weak, since it does not impose constraints on the behavior of the covariance matrix Sx = ^X;n 

428 in the sequence of distributions Pn (tt, G N), except for strict positive definiteness for all n. Con- 

429 dition (D5,(ii)) excludes perfect collinearity, where the fixed upper bound on partial correlations 

430 places some additional restrictions on the design. Regarding the conditional distribution of Y 

431 given X, we require partial faithfulness. This becomes more expUcit by invoking Theorem 1: 

432 partial faithfulness follows by assuming condition (C2) in Section 2 for every n, which involves 
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Example 4. Consider a Gaussian linear model as in (1) with p = 4, pefF = 3, cr^ = 1, jjlx 

(o,...,of 



pi = -0.4, p2 = 0.2, 



433 the regression coefficients only. Conditions (D4) and (D5,(i)) place additional restrictions on 

434 both the design X and the conditional distribution of Y given X. 

435 Assumption (D4) is used for controlling the type 11 errors in the many tests of the PC-simple 

436 algorithm, see the proof of Theorem 4. This assumption is slightly stronger than requiring that all 

437 non-zero regression coefficients are larger than a detectability-threshold, which has been previ- 

438 ously used for analyzing the Lasso in Meinshausen & Biihlmann (2006), Zhao & Yu (2006) and 

439 Meinshausen & Yu (2009). Clearly, assumptions on the design X are not sufficient for consistent 

440 variable selection with any method and some additional detectability assumption is needed. Our 

441 assumption (D4) is restrictive, as it does not allow small non-zero low-order partial correlations. 

442 Near partial faithfulness (Robins et al, 2003), where small partial correlations would imply that 

443 corresponding regression coefficients are small, would be a more realistic framework in practice. 

444 However, this would make the theoretical arguments much more involved, and we do not pursue 

445 this in this paper. 

446 Although our assumptions are not directly comparable to the neighborhood stabihty or irrep- 

447 resentable condition for the Lasso, it is easy to construct examples where the Lasso fails to be 

448 consistent while the PC-simple algorithm recovers the true set of variables, as shown by the 

449 following example. 
450 
451 
452 
453 
454 
455 
456 
457 

458 where Pi, P2, /?3 are fixed i.i.d. realizations from A/'(0, 1) and P4 = 0. 

It is shown in Zou (2006, Cor. 1) that the Lasso is inconsistent for variable selection in this 

460 model. On the other hand, (Dl) holds because of Theorem 1, and also (D5) is true. Since the 
dimension p is fixed, (D2), (D3) and (D4) hold automatically. Hence, the PC-simple algorithm 
is consistent for variable selection. It should be noted though that the adaptive Lasso is also 
consistent for this example. 

464 

465 We can slightly modify Example 4 to make it high-dimensional. Consider pefF = 3 active vari- 

466 ables, with design and coefficients as in Example 4. Moreover, consider p„ — pcff noise co- 

467 variates which are independent from the active variables, with p^ satisfying (D2). Let the design 

468 satisfy (CI) and (D5,(ii)), for example by taking the noise covariates to be mutually independent. 

469 Then assumptions (D1)-(D5) hold, implying consistency of the PC-simple algorithm, while the 

470 Lasso is inconsistent. 
471 

4-72 5 '3- Asymptotic behavior of correlation screening 

473 Correlation screening is equivalent to sure independence screening of Fan & Lv (2008), but our 

474 assumptions and reasoning via partial faithfulness are very different from the work of Fan & Lv. 

475 Denote by An^ {a) the correlation screening active set, estimated from data, using significance 

476 level a, obtained from the first step of the sample version of the PC-simple algorithm. We do not 

477 require any sparsity conditions for consistency. We define: 
478 

479 (El) as assumption (D4) but for marginal correlations cor(y, X'^^') = pn{Y,j) only. 

430 (E2) as assumption (D5) but for marginal correlations cor(y, X^^) = pn{Y,j) only. 



( I pi Pi Pi\ 

Pi 1 Pi P2 
Pi Pi 1 P2 
\P2 P2 P2 1 / 



Partial faithfulness and the PC-simple algorithm 



11 



pr0^\an) D A} = 1 - 0{eM-Cn^''"^)} ^ 1 (n ^ oo), 



Numerical results 



481 Theorem 5. Consider the linear model (1) and assume (Dl), (D2), (El) and (E2). Then 

482 there exists a sequence q;„ ^ (n — oo) and a constant C > such that: 
483 
484 

485 where d > is as in (El). 

486 ^ ,„ 

4^-7 A proof is given in the Appendix. A possible choice for q;„ is a„ = 2{1 — ^{n ' Cn/2)}. As 

4gg pointed out above, we do not make any assumptions on sparsity. However, for non-sparse prob- 

439 lems, many correlations may be non-zero and ^1^1 can still be large, for example almost as large 

49Q as the full set {1, . . . 

492 Under some restrictive conditions on the covariance Sx of the random design. Fan & Lv 

492 (2008) have shown that correlation screening, or sure independence screening, is overestimating 

493 the active set A, as stated in Theorem 5. Theorem 5 shows that this result also holds under very 

494 different assumptions on T,x when partial faithfulness is assumed in addition. Hence, our result 

495 justifies correlation screening as a more general tool than what it appears to be from the setting 
495 of Fan & Lv (2008), thereby extending the range of applications. 

497 
498 
499 

500 6 1. Analysis for simulated data 

501 We simulate data according to a Gaussian Unear model as in (1) with 6 = and p covariates 

502 with nx = {0, . . . , 0)^ and covariance matrix T.x-i,j = p'*^-'', where denotes the (i,i)th 

503 entry of S. In order to generate values for /?, we follow (C2): a certain number pcff of coefficients 

504 Pj have a value different from zero. The values of the nonzero /3jS are sampled independently 

505 from a standard normal distribution and the indices of the nonzero Pjs are evenly spaced between 

506 1 and p. We consider two settings: 
507 
508 
509 

510 We evaluate the performance of the methods using receiver operating characteristic curves 

511 which measure the accuracy for variable selection independently from the issue of choosing 

512 good tuning parameters. We compare the PC-simple algorithm to the Lasso (Efron et al., 2004) 

513 and Elastic Net (Zou & Hastie, 2005), using the R-packages pcalg, lars and elasticnet, 

514 respectively. For Elastic Net, we vary the £^-penalty parameter only while keeping the £^ -penalty 

515 parameter fixed at the default value from the R-package. 

516 In the low-dimensional settings shown in Figures 2(a), 2(c), 2(e), the PC-simple algorithm 

517 clearly dominates the Lasso and Elastic Net for small false positive rates, which are desirable 

518 in many applications. When focusing on the false positive rate arising from the default value 

519 for a=0-05 in the PC-simple algorithm, indicated by the vertical lines, the PC-simple algorithm 

520 outperforms the Lasso and Elastic Net by a large margin. If the correlation among the covariates 

521 increases, the performance of Elastic Net deteriorates, whereas the performances of the PC- 

522 simple algorithm and the Lasso do not vary much. 

523 In the high-dimensional settings shown in Figures 2(b), 2(d), 2(f), the difference between the 

524 methods is small for small false positive rates. The Lasso performs best. Elastic Net is worst, 

525 and the PC-simple algorithm is somewhere in between. For larger false positive rates, the differ- 

526 ences become more pronounced. Up to the false positive rate corresponding to the default value 

527 of q;=0 05, the PC-simple algorithm is never significantly outperformed by either the Lasso or 

528 Elastic Net. 



low-dimensional: p = 19, peff = 3, n = 100; p e {0,0-3,0-6} with 1000 replicates 
high-dimensional: p = 499, pefF = 10, n = 100; p G {0,0-3,0-6} with 300 replicates 
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We refer to the working paper "Stability selection" by N. Meinshausen and P. Biihlmann for a 

more principled way to choose the tuning parameter a. Further examples, with p = 1000, pcff = 
5, n = 50 and equi-correlated design ^X;i,j = 0-5 for i ^ j and = 1 for all i, wee reported 

in Biihlmann (2008). 

The computing time of the PC-simple algorithm on 10 different values of a has about the 
same order of magnitude as the Lasso or Elastic Net for their whole solution paths. Hence, the 
PC-simple algorithm is certainly feasible for high-dimensional problems. 



OM 



t:b5 otw ots 

False positive rate 
(a) Low dimensional, p = 0. 




0.000 0.002 



0.004 0506 0.008 
False positive rate 

(b) High dimensional, p = 0. 




530 



OTIO 0715 
False positive rate 

(c) Low dimensional, p=0-3. 



T50 



0^000 0.002 



"0304 0.006 
False positive rate 



OTO 5310 



(d) High dimensional, p=0-3. 




530 



T35 OTIO 5715" 
False positive rate 

(e) Low dimensional, p=0-6. 



050 



0300 5302" 



0.004 356 5358 5310 
False positive rate 



(f) High dimensional, p=0-6. 



Fig. 2. Receiver operating characteristic curves for the sim- 
ulation study in Section 6-1. The solid line corresponds to 
the PC-simple algorithm, the dashed line to the Lasso and 
the dotted line to Elastic Net. The solid vertical lines indi- 
cate the performance of the PC-simple algorithm using the 
default a=0-05. 



6-2. Prediction optimal tuned methods for simulated data 
We now compare the PC-simple algorithm to several existing methods when using predic- 
tion optimal tuning. It is known that the prediction-optimal tuned Lasso overestimates the true 
model (Meinshausen & Biihlmann, 2006). The adaptive Lasso (Zou, 2006) and the relaxed Lasso 
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577 (Meinshausen, 2007) correct Lasso's overestimating behavior and prediction-optimal tuning for 

578 these methods yields a good amount of regularization for variable selection. 

579 We simulate from a Gaussian linear model as in (1) with p = 1000, peff = 20,n = 100, and 

If^ 5 = 0, MX = (0, . . . ,0)^, = 0-5l^-^l, = 1, 
582 



. . . , /320 i-i-d. ~ A/'(0, 1), /32i = . . . = Aooo = 0, 



(11) 



argnun^gRp 



J2{Y,-Xrpf + Xj2wj'\P,\ 

1=1 7=1 



583 using 100 replicates. We consider the following performance measures: 
584 

585 II/? - PWl = E^=i(/3, - (^jf (MSE Coeff) 

586 Ex[{XT0 - = 0- /3)cov(X)(/3 - f5f (MSE Pred) 

587 J2%i 103 + 0, /3, + 0)/ Y]]=i 103 + 0) (true positive rate) 

588 y:%i ^ 0' = 0)/ ELi 103 = 0) (false positive rate) 
589 

590 where /(•) denotes the indicator function. 

59 \ We apply the PC-simple algorithm for variable selection and then use the Lasso or the adaptive 

592 Lasso to estimate the coefficients for the sub-model selected by the PC-simple algorithm. We 

593 compare this procedure to the Lasso, the adaptive Lasso and the relaxed Lasso. For simplicity, 

594 we do not show results for Elastic Net, since this method was found to be worse in terms of 

595 receiver operating characteristic curves than the Lasso, see Section 6-1. 

595 Prediction optimal tuning is pursued with a validation set having the same size as the training 

597 data. For the adaptive Lasso, we first compute a prediction-optimal Lasso as initial estimator 

59g Anit, and the adaptive Lasso is then computed by solving the following optimization problem: 

599 
600 
601 
602 

603 where wj^ = \(3init,j\~'^ and A is again chosen in a prediction-optimal way. The computations 

604 are done with the R-package lars, using re-scaled covariates for the adaptive step. The relaxed 

605 Lasso is computed with the R-package r e 1 ax o . The PC-simple algorithm with the Lasso for es- 

606 timating coefficients is computed using the R-packages pcalg and lars, using optimal tuning 

607 with respect to the a-parameter for the PC-simple algorithm and the penalty parameter for Lasso. 

608 For the PC-simple algorithm with the adaptive Lasso, we first compute weights Wj as follows. 

609 If the variable has not been selected, we set Wj = 0. If the variable has been selected, we let Wj 

610 be the minimum value of the test statistic (n — 3 — \S\Y/'^Ziy,j \ S) over all iterations of the 

611 PC-simple algorithm. With these weights Wj, we then compute the adaptive Lasso as defined 

612 above. 

613 The results are shown in Figure 3. As expected, the Lasso yields too many false positives, 

614 while the adaptive Lasso and the relaxed Lasso have much better variable selection properties. 

615 The PC-simple based methods clearly have the lowest false positive rates, but pay a price in 

616 terms of the true positive rate and mean squared errors. In many applications, a low false positive 

617 rate is highly desirable even when paying a price in terms of power. For example, in molecular 

618 biology where a covariate represents a gene, only a limited number of selected genes can be 

619 experimentally vahdated. Hence, methods with a low false positive rate are preferred, in the 

620 hope that most of the top-selected genes are relevant, as sketched in the next section. 
621 

622 6-3. Real data: riboflavin production by Bacillus subtilis 

623 We consider a high-dimensional real data set about vitamin riboflavin production by the bac- 

624 terium B. subtilis, kindly provided by DSM Nutritional Products. There is a continuous re- 
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634 
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638 |2- . . . |oJ : _ 

639 ^ ^ 

640 " 
641 
642 

643 Fig. 3. Boxplots of performance measures for the simula- 

544 tion study in Section 6-2 considering the following predic- 

tion optimal tuned methods: the PC-simple algorithm with 
Lasso coefficient estimation (pel), the PC-simple algorithm 

646 with adaptive Lasso (peal), the adaptive Lasso (al), the re- 

647 laxed Lasso (r) and the Lasso (1). 

648 
649 

sponse variable Y which measures the logarithm of the riboflavin production rate, and there 
J aiep = 4088 covariates corresponding to the logarithms of expression levels of genes. The main 

go^l is to genetically modify B. subtilis in order to increase its riboflavin production rate. An 
important step to achieve this goal is to find genes which are most relevant for the production 
654 ""^te. 

We use data from a genetically homogeneous group of n = 71 individuals. We run the PC- 
simple algorithm on the full data set for various values of a. Next, we compute the Lasso and 
Elastic Net, choosing the tuning parameters such that they select the same number of variables 
as the PC-simple algorithm. 

Table 1 shows that there is overlap between the selected variables of the three different meth- 
ods. This overlap is highly significant when calibrating with a null-distribution which consists of 
random noise. On the other hand, we see that the variable selection results of the Lasso and Elas- 
ti'^ Net are more similar than the results of the PC-simple algorithm and either of these methods. 
Hence, the PC-simple algorithm seems to select genes in a different way than the penalty-based 
methods Lasso and Elastic Net. This is a desirable finding, since for any large-scale problem, 

665 want to see different aspects of the problem by using different methods. Ideally, results from 

different methods can then be combined to obtain results that are better than what is achievable 

gg-y by a single procedure. 

668 

559 Appendix 

670 Proofs 

671 Proof of Theorem 1. Consider the linear model (1) satisfying assumptions (CI) and (C2). In order to 

672 prove that the partial faithfulness assumption holds almost surely, it suffices to show that the following 
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673 
674 
675 
676 
677 
678 
679 
680 
681 
682 
683 
684 
685 
686 
687 
688 
689 
690 
691 
692 
693 
694 
695 
696 
697 
698 
699 
700 
701 
702 
703 
704 
705 
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707 
708 
709 
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711 
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713 
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715 
716 
717 
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720 



a for PC-simple 


Selected 


PC-Lasso 


PC-Enet 


Lasso-Enet 


0001 


3 








2 


001 


4 


2 


1 


3 


0-05 


5 


2 


1 


3 


015 


6 


3 


2 


3 



Table 1. Variable selection for a real data set on riboflavin production by B. subtilis. The 
columns 2 to 5 show the number of variables selected by the PC-simple algorithm, the number 
of variables selected by both the PC-simple algorithm and the Lasso, the number of variables 
selected by both the PC-simple algorithm and Elastic Net, and the number of variables that were 
selected by both the Lasso and Elastic Net. 



holds for all j G {1, . . . ,p)anAS Q {j}'^: /3j 7^ impUesthat X^-?) | X^-^)) 7.^ almost surely with 
respect to the distribution generating the /3jS. 

Thus, let j e {1, ... ,39} such that Pj ^ 0, and let 5* C {j}<^. We recall that p{Y, | X^^^) = if 
and only if the partial covariance parcov(Y, X^^^ \ X^^^^) between Y and X^^^ given X^^^ equals zero 
as can be seen in Anderson (1984, page 37, definition 2-5-2). Partial covariances can be computed using 
the recursive formula given in Anderson (1984, page 43, equation (26)). This formula shows that the 
partial covariance is Unear in its arguments, and that parcov(e, X^^^ | X^^^) = for all j G {1,. . . ,p} 
and,S C {j}<^. Hence, 



parcov(y, X(^) | X'-^^) = parcov((5 + -\- e,X<^^'> \ X^^^) 

r=l 

P 

= ^/3^parcov(XW,X(j) | X^^^) 

r=l 

P 

= /3jparcov(X(^\x(j') | X^^^) + ^ ^,parcov(X(''\ X^^'^ | X^^^). 

r—l.r^j 

Since /3j ^ by assumption, and since parcov(X(-'^ , X^^^ \ X^'^^) ^ by assumption (CI), the only way 
for parcov(y, X^^^ | X^^^ ) to equal zero is if there is a special parameter configuration of the /J^s, such 
that 



^ /3^parcov(x(''),x(j') | x'-^^) = -PjpavcoY{X^^\ X'-J^ \ X^^^). 
But such a parameter constellation has Lebesgue measure zero under assumption (C2). 



(Al) 
□ 



Proof of Corollary 1. The implication from the left to the right hand side follows from the fact that 

/3j 7^ in the Hnear model (1) if and only if p{Y, X^^^ \ X({J>'^)) ^ 0. The other direction follows from 
the definition of partial faithfulness, by taking the negative of expression (2). □ 

Proof of Theorem 2. Suppose that {X, Y) = {X'^^\ . . . , X^p\Y) is linearly F-faithful to a directed 
acyclic graph G in which Y is childless, i.e., any edges between Y and the X^-'^s, j = 1, . . . ,p, point 
towards Y . We will show that this implies that the distribution of (X, Y) is partially faithful, by showing 
that p(r, X(J) I X«J>^)) / implies that p{Y, X^^^ \ X'^^^) ^ for all S C {j}"^. 

Thus, let j e {1, . . . ,p} such that p{Y,X^i'^ \ X'^^^^' ^) / 0. By linear F-faithfulness, this implies that 
Y and X^^) are not d-separated by X(^J> ) in G, meaning that X'^^^^ does not block all d-connecting 
paths between X^^'^ and Y. All paths between X^-'^ and Y must be of the form X^^^ — — 
— > Y, where — denotes an edge of the form <— or First suppose that r ^ j. Then, because X^''^ 
cannot be a colUder on the given path (since we know that the edge from X^''^ to Y points towards Y), the 
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721 path is blocked by X^'') e X^^^'i"\ and hence the path is blocked by X(^J>'^). Thus, since X({-'>°) does 

722 not block all paths between X^^^ and Y, there must be a path where r = j, or in other words, there must 

723 be an edge between X'^^^ and Y: X^^^ — > Y. Such a path X'^^^ — » Y cannot be blocked by any set 

724 S C {j}'" . Hence, there does not exist a set <S that d-separates X^^^ and Y. By Unear y-faithfulness, this 

725 impUes that p{X'^^\Y | X^-^)) ^ for all 5 C {j}^. □ 

Proof of Theorem 3. By partial faithfulness and equation (8), A C ^['"reach] Hence, we only need to 

"^■^"^ show that A is not a strict subset of ^[™roach] -y^e do this using contra-position. Thus, suppose that 

728 j[ (2 ^[™roach] strictly. Then there exists a j G _4,['""'°^°''' such that j ^ yl. Fix such an index j. Since 

729 j g ^[mreaoh]^ know that 
730 

731 p(y, ^^^^ I X^^^) ^ for all S C ^[™-ch-i] \ {j-} ^jth |5| < .^^^^^^^ _ i. (A2) 

732 xijjs statement for sets S with |5| = mj-cach — 1 follows from the definition of iteration mrcach of the 

733 PC-simple algorithm. Sets S with lower cardinality are considered in previous iterations of the algorithm, 

734 and since At^l D D . . ., all subsets S C ^[™rcaoh-i] with \S\ < TOreach - 1 are considered. 

735 We now show that we can take S = Ain (A2). First, the supposition A C ^['""aoh] our choice of 

736 j imply that 



763 
764 



737 
738 

739 Moreover, A C ^[™'each] impUes that |^| < |>l[™reach]| _ i Combining this with |^["'"ach]| < mreach 

740 yields that |^| < mreach — 1- Hence, we can indeed take <S = .A in (A2), yielding that p{Y,X^^'' \ 

741 ^ 0. 

742 On the other hand, j ^ A implies that Pj = 0, and hence p{Y, X^^^ \ X^"^^) = 0. This is a contradic- 

743 tion, and hence A cannot be a strict subset of ^i^wach] □ 

744 Proof of Theorem 4. A first main step is to show that the population version of the PC-simple algo- 

745 rithm infers the true underlying active set An, assuming partial faithfulness. We formulated this step as a 

746 separate result in Theorem 3 . 

747 The arguments for controlling the estimation error due to a finite sample size are similar to the ones 

748 used in the proof of Theorem 1 in Kalisch & Biihlmann (2007). We proceed in two steps, analyzing first 

749 partial correlations and then the PC-simple algorithm. 
750 

We show an exponential inequality for estimating partial correlations up to order m„ = o{n). We use 
the following notation: i^Tj"" = {5 C {0, . . . ,p„} \ {j}; \S\ < nin} = 1,... ,Pn)- We require more 

752 general versions of assumptions (D4) and (D5) where the cardinaUty of the condition sets are bounded by 

753 the number to„ : 
754 

755 (D4„„) The partial correlations Pn{Y,j \ S) = | X^'^)) satisfy: 

756 . 

757 inf ||pn(^,i I 5)1; j = 1,... ,Pn, S C {j}'", \S\ < nin with pn{Y,j | 5) 7^ j > c„, 

758 . 

rjcq where c„ = 0(n'') for some < d < b/2, and 6 is as in (D3). 

(D5„„) The partial correlations p„{Y,j \ S) and p„(i, j | S) = p(X('), X^^) | X^'^)) satisfy: 

761 (i) sup \pn{Y,j\S)\<M <1, {ii) sup \pn{i,j\S)\<M <1. 

762 nj,5C{j}C,|5|<m„ ",i7^i,<SC{iJ}C,|S|<m„ 



We will later see in Lemma 3 that we need m„ < peff „ only, and hence, assumptions (04^^) and (D5m,„) 
coincide with (D4) and (D5), respectively. 
765 We have, for m„ < n - 4 and < 7 < 2, 



766 
767 

768 seKr\3=i,...,Pn V4 + 72 



'^^'^ sup pr{ |p„ (Y, j\S)- pn (r, J I 5) I > 7} < Cm exp(n - m„ - 4) log ( 
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769 where < Ci < oo depends on M in (D5to„) only. This bound appears in Kalisch & Buhlmann (2007, 

770 Corollary 1): for proving it, we require the Gaussian assumption for the distribution and (D5m„)- It 

77 1 is now straightforward to derive an exponential inequality for the estimated Z-transformed partial cor- 

772 relations. We define Zn{Y,j \ S) = g{pn{yj I ^)} Zn{Y,j \ S) = g{pn{Y,j \ S)}, where g{p) = 

773 ilog{(l + p)/(l-p)}. 

Lemma 1 . Suppose that the Gaussian assumption from (Dl) and condition (D5„i^) hold. Define L = 
775 1/{1 - (1 + M)2/4}, with M as in assumption (D5m, )■ Then, for mn < n - A and < 7 < 2L, 

776 

777 sup pr{\Zn{Y,j \ S) - Zn{Y,j \ S)\ > 7} 



778 



794 
795 



810 
811 



•SGifT"" ,j=l,...,p, 

•4-(7/i)^ 



(n - 4 - rUn) log 



4+(7/i)2 



-exp{-C2(n-m„)} 



779 <0{n){ exp 

780 V 

78 1 for some constant C2 > 0. 
782 

783 We omit the proof since this is Lemma 3 in Kalisch & Buhlmann (2007). 

784 

We now consider a version of the PC-simple algorithm that stops after a fixed number of m iterations. 

78^ If m > mreach, where TOrcach is the estimation analogue of (9), we set ^I™! = ^["*reach]_ -s^^ denote this 

786 version by PC-simple(m) and the resulting estimate by A{a, m). 

787 

788 Lemma 2. Assume (D1)-(D3), (D4m„) and (D5„i^). Then, for m„ satisfying nin > T7ireach,n ctnd 

"789 TTT-n = 0{n^~^) with b as in (D3), there exists a sequence such that 

790 
791 

792 A concrete choice of a„ is a„ = 2{1 — $(n^/^c„/2)}, where c„ is the lower bound from (04^^), which 

793 is typically unknown. 



pr{An{oin,'nin) = An} = 1 - 0{exp(— Cn )} — *■ 1 (n — > oo) for some C > 0. 



<A U ^.15 <0(p™"+i) sup vAE^\s), (A3) 



Proof. Obviously, the population version of the PC-simple(TO„) algorithm is correct for m„ > 
mreach,n, scc Theorem 3. An error can occur in the PC-simple(m„) algorithm if there exists a covari- 
ate X'^^^ and a conditioning set S G ifj"" for which an error event Sj|5 occurs, where Ej\s denotes the 
7^7 event that an error occurred when testing pn{Y. j \S) = 0. Thus, 

798 

799 pr{an error occurs in the PC-simple(TO„)-algorithm} 

800 
801 
802 

803 using that the cardinaUty of the index set {<S G K™",j = 1, . . . in the union is bounded by 

804 0(j3;^"+i).Now 
805 

806 E^\s = El^s^El{s, (A4) 

807 
808 

809 type I error : (n - \S\ - 3)^/^\Zn{Y,j \ S)\ > - a/2) and Zn{Y,j \ S) = 0, 

type n error E'.(s : (n - \S\ - 3)^/^\Zn{Y,j \ S)\ < - a/2) and Zn{Y,j \S)j^O. 



where 



812 Choose a = a„ = 2{1 - <i>{n^/^Cn/2)}, where c„ is from (D4„ ). Then, 
813 

814 sup pr(4|5)= sup pr\\Zn{Y,j\S)-Zn{Y,j\S)\>{n/{n-\S\-3)y/^Cn/2 

815 5eK7"',j=i,...,p„ seK^",j >- 

816 < 0(n)exp{-C3(n-m„)c2}, (A5) 
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817 for some C3 > 0, using Lemma 1 and the factthat log{(4 - 52)/(4 + ^2-)| < -j2/2for0 < 6 < 2.Fur- 

818 thermore, with the choice of a = above, 

819 r 

020 sup V<E]{s) = sup pr |Z„(y,i I 5)1 < {n/(n - \S\ - 3)}V2c„/2 

5e/f'"",j = l,...,p„ SGK"^",j 

821 

822 < sup pr(\Zr,{Y,j\S)-Zn{Y,j\S)\>Cn[l-{n/{n-\S\-3)y/y2]), 

823 ^^^^""'■^^ 

824 because inf5e;^™„,^{|z„(r, j | S)\; Zn{Y,j | 5) ^ 0} > c„ since \g{p)\ = |i log{(l + p)/{l - p)}\ > 
§25 IpI for aU p and using assumption {D4m„). This shows the crucial role of assumption (04^^) in controlling 
826 the type II error. By invoking Lemma 1 we then obtain: 

sup pt{E¥s) < Oin) exp{-C4(n - m„)4} (A6) 



827 

828 56K""/'^"'' 
829 

g3Q for some C4 > 0. Now, by (A3)-(A6) we get 

§^ ^ pr{an error occurs in the PC-simple(m„)-algorithm} 

< 0\p^-+'n exp{-C5(n - m„)c2 }] < 0[n-('"»+i)+i exp{-C5(n - m„)n-2'^}] 
334 = O [exp {a(m„ + 1) log(n) + log(n) - Csin^'^'^ - m„n-2'^)}] = o(l), 

835 because n^~'^'^ dominates all other terms in the argument of the exp-function, due to m„ = 0{n^~^) and 

836 the assumption in (1)4^^) that d < 6/2. This completes the proof. □ 
837 

838 Lemma 2 leaves some flexibiUty for choosing m„. The PC-algorithm yields a data-dependent stopping 

839 level mreach,n> that is, the sample version of (9). 
840 
841 

842 Pr{mreach,n = ^Tlreach.n) = 1 - 0{exp(-Cn^~^'^)} ^ 1 (n ^ Oo) 

for some C > 0, with d is as in (D4). 

844 

g45 Proof. Consider the population version of the PC-simple algorithm, with stopping level mreach as de- 

846 fined in (9). Note that TOicach = 'Tiroach.n = 0(n^~^) under assumption (D3). The sample PC-simple(TO„) 

algorithm with stopping level in the range of m„ > Wreach (to„ = 0(n^^'')), coincides with the popu- 
lation version on a set A having probabiUty P[A] = 1 — 0{exp(— Cn^"^"^)}, see the last formula in the 
proof of Lennma 2. Hence, on the set A, mreach.n = 'Tireach- D 



Lemma 3. Assume (D1)-(D5). Then, 



848 
849 

850 Lemma 2 with m„ = peff„ together with Lemma 3, and using that rrireach.n < peffn' complete the 

85 1 proof of Theorem 4. □ 

Proof of Theorem 5. By definition. An C ^'^J for the population version. Denote by Zn{Y,j) the 

853 quantity as in (10) with 5 = and by Zn{Y,j) its population analogue, i.e., the Z-transformed popu- 

854 lation correlation. An error occurs when screening the jth variable if Z„ (Y, j) has been tested to be zero 

855 but in fact z„(F, j) ^ 0. We denote such an error event by Ej^ . Note that 

856 „ 

857 sup pi{Ej ) <0{n)exp{-CincJ, 

j=i,...,p„ 

858 

859 some Ci > 0, see formula (A6) above. We do not use any sparsity assumption for this derivation, but 

860 invoke (El) which requires a lower bound on non-zero marginal correlations. Thus, the probability 
of an error occurring in the correlation screening procedure is bounded: for some C2 > 0, 

861 

862 pr( Uj=i,...,p„ £;f ) = 0(p„n) cxp{~Cincl) = 0[exp{(l + a) log(n) - Cm^-^"^}] 

863 = 0{eM-C2n^~^'^)}. □ 
864 
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